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A set of direct simulations of zero-pressure gradient, turbulent boundary layer flows 
are conducted using various span widths (62 — 630 wall units), to document their influence 
on the generated turbulence. The FDL3DI code that solves compressible Navier-Stokes 
equations using high-order compact-difference scheme and filter, with the standard re- 
cycling/rescaling method of turbulence generation, is used. Results are analyzed at two 
different Re^ values (500 and 1,400), and compared with spectral DNS data. They show 
that a minimum span width is required for the mere initiation of numerical turbulence. 
Narrower domains (< 100 w.u.) result in relaminarization. Wider spans (> 600 w.u.) are 
required for the turbulent statistics to match reference DNS. The upper-wall boundary 
condition for this setup spawns marginal deviations in the mean velocity and Reynolds 
stress profiles, particularly in the buffer region. 

I. Introduction 

Turbulent boundary layers have been the subject of intensive research for many decades, partly be- 
cause of their pervasive nature. With advances in computing resources, extending the state-of-the-art in 
understanding these complex flows has been further enabled. 

Ever since Spalart’s seminal paper, in which he presented DNS of turbulent boundary layer on a flat plate 
for a range of Re# values (225 — 1, 410), successive efforts have managed to extend this high Reynolds number 
limit. Rai and Moin used a high-order finite- difference approach to study transition and turbulence of a 
spatially evolving boundary layer up to Re# = 2, 300. Ferrante and Elghobashi used the recycling/rescaling 
technique to study spatially developing incompressible boundary layers up to Re# = 1,430. 

More recently, a number of research efforts have made significant advances in extending that limit further. 
Lee and Sung 4 also used a similar technique as this paper, and computed turbulent boundary layer up to 
Re^ = 2,500. Schlatter and Orlii were able to further extend DNS data to Re# = 4,300 using spectral 
methods. Within the last few years, Sillero et al. managed to breach the 6,000 limit and generate results 
up to Re <9 = 6, 650. Borrell et al. furthered this by computing up to Re^ = 6, 800. 

All of these efforts have been restricted to the incompressible regime, and many of them used spectral 
methods. The applicability of such techniques to practical flow situations is rather limited. Not necessarily 
unique, this study steps outside that usual realm by using a finite-difference code that solves compressible 
Navier-Stokes equations. Documenting challenges encountered, and subsequent solutions, will be useful 
for future simulations of compressible turbulent flows. This study reports the requirements for accurate 
predictions of turbulent boundary layers, including the minimum spanwise width of the domain. In that 
sense, this work follows Jimenez and Moin’s article on the impact of width on turbulent channel flow. 
In their study, they conducted DNS of low Reynolds number channel flow to study the morphology and 
dynamics of the most basic flow unit of turbulence. The fundamental period of intermittency was identified, 
and the impact of different domain widths on the growth of turbulence was investigated for the unsteady 
channel flow problem. This study analyzes turbulent boundary layer on a flat plate instead. 

This study is a part of an overall effort to understand turbulent, separated flows. There continues 
to be a dearth in the physics knowledge-base for such phenomena. This is reflected in the computations of 
separated flows using RANS models, which significantly deviate from experimental observations, particularly 
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in reattachment and recovery. The understanding that this study presents will aid to further improve the 
state-of-the-art in flow predictions using relatively inexpensive RANS and hybrid RANS/LES methods. 
As a precursory exercise, it is important to establish the capacity to predict attached turbulent boundary 
layers at various Reynolds numbers. In the process, requirements for accurate predictions of the physics are 
parameterized. Included amongst these are the minimum width of span and the ideal mesh resolution at 
the wall. Results from this effort are being used in a simultaneous study of turbulent, separated flows past 
a backward- facing step and a wall-mounted hump. They are also being used to set up simulations at higher 
Re# values. 

The paper is structured as follows: Section II briefly explains the numerical methodology used for the 
simulations. Section III documents the results for both Re# = 500 and 1,400 simulations. Conclusions are 
drawn in Sec. IV. 


II. Numerical Method 

II. A. Governing Equations, Compact Difference Scheme and Filter 

The finite difference code (a variant of AFRL’s FDL3DI 5 ) used in this study solves the compressible Navier- 

Stokes equations for an ideal gas. The sixth-order compact-difference scheme of Lele 1 is used to solve the 
governing equations in transformed curvilinear coordinates. 2,13 Fourth-order explicit Runge-Kutta (RK-4) 
is used for time integration. During each time step, an eighth-order low-pass spatial filtering scheme is 
applied to the conservative variables to ensure stability, along with second- and fourth-order near-boundary 
formulations. A filtering optimization parameter of otf = 0.495 is set. To ensure stability and accuracy 
of time integration, a Courant-Friedrichs-Lewy (CFL) number in the vicinity of 0.4, based on the acoustic 
speeds, is used in the RK-4 scheme. 

This code was previously designed to solve for supersonic flows, specifically those involving shock- wave 
and boundary layer interaction problems. Therefore, the localized artificial diffusivity (LAD) method of 
shock-capturing was previously installed by Kawai et ah, ’ and was used in a detailed study of the flow 
physics of such interactions by Morgan et al. This method adds artificial coefficients to the fluid transport 
terms, and takes effect in the vicinity of shocks. Since the flow regime under consideration is subsonic, LAD 
is not included in the current study. Thus, other than the low-pass filter, no additional sub-grid scale model 
is used. In that sense, this effort can be considered implicit-LES, or more precisely “under-resolved DNS.” 

II. B. Recycling/Rescaling Turbulence Inflow Method 

For generating the turbulent inflow, the standard recycling/rescaling method of Lund et al. is used in 
the code. Urbin and Knight 1 extended this method to compressible flows for solving a Mach 3 turbulent 
boundary layer. This was implemented as a part of the effort to study the interaction of an oblique shock 
wave and a turbulent boundary layer by Morgan et al. The approach extracts the instantaneous velocity 
profile from a plane downstream of the inflow, rescales and reintroduces it as the inflow. The schematic 
of the procedure is shown in Fig. 1. It has been shown to produce a realistic turbulent boundary layer, 
yielding statistics that are in good agreement with experiment and DNS. To further prevent contamination 
of the solution with spurious spatio-temporal correlations generated by this reintroduction procedure, a non- 
constant reflection of the recycled turbulence plane is applied at randomly distributed time intervals. This 
method, called recycling/rescaling with dynamic reflection (RR+DR), has been suitably adjusted to solve 
subsonic compressible flows. 
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III. Results 


The setups for validating the code for subsonic turbulent boundary layers are detailed here. A separate 
simulation was required for each Re#. 

For both simulations, a domain of streamwise length 15 Jr and height 3 Jr was used, where Jr is the 
desired boundary layer thickness at the velocity capture plane. This capture plane was located about 12 Jr 
downstream of the inlet (Fig. 1), and was also the location of the target Re#. The initial mesh in each 
case had a spanwise width of |Jr, and this was doubled successively as required. This mesh consisted of 
500 x 100 x 27 points, and was uniformly spaced along the streamwise and spanwise directions. A hyperbolic 
tangent stretching was set in the wall normal direction. 

For the boundary conditions, the inlet was computed during the recycling/rescaling procedure, and a 
freestream pressure outflow was set for the outlet. A no-slip, adiabatic wall was used for the bottom wall. 
First-order extrapolation of the state variables was performed at the top wall. Along the spanwise direction, 
periodic boundary conditions were used. 

For both cases, sea level conditions were assumed, with a freestream Mach number of 0.3. A 0.1% random 
perturbation was included as a part of the initial condition. 

III. A. Re# = 500 

To simulate Re# = 500 at the location 12 Jr downstream of the inlet, Re# at the inlet was set to 350, based 
on Prandtl’s one-seventh power law approximation . 

At these flow conditions, the mesh resolution at the wall was computed to be Ax + = 7.5 and Az + = 2.3, 
along the streamwise and spanwise directions, respectively. The overall spanwise width of the domain was 
L z = 62 w.u. In the wall normal direction, the spacing was A y + = 0.25 at the wall, and was stretched using 
a hyperbolic-tangent function. With such a resolution, this mesh is almost twice as fine as that used by 
Spalart for his spectral DNS. But, as Moin and Mahesh noted, 22 spatial resolution requirements are dictated 
by the choice of numerical schemes. They showed that sixth-order Pade requires 60% finer meshes than 
those demanded in Fourier spectral methods. 

Flow development was monitored in terms of u velocity component contours, to determine the initial- 
ization of turbulence. As shown in the time series in Fig. 2, after approximately 12 flow through time 
cycles (FTTs), the laminar nature of the boundary gradually starts to break down. After about 15 FTTs, 
wavy Kelvin-Helmholtz-like instabilities are evident in the domain. These disturbances, by means of the 
recycling/rescaling process, are fed back into the domain, resulting in an increase in their intensity. This 
process continues to perturb the flow until about 22 FTTs, when the flow starts to go turbulent. At this 
point, instead of sustaining the turbulent motion in the domain, the flow relaminarizes and the disturbances 
get dampened out. It should be noted that this phenomenon is a purely numerical transition because of the 
recycling process, and does not correspond to the physical transition from a laminar to a turbulent profile. 

A similar observation was made by Jimenez and Moin for channel flows, where domains with spanwise 
widths lesser than a threshold failed to sustain turbulent motion. In their estimate, for spans narrower 
than 100 w.u., laminar solutions alone were found. This observation is in agreement with their estimate. 
Periodicity on such narrow domains compromises 2-point correlations, damping out the initiated turbulence. 
To confirm this, the L z = 62 w.u. domain was doubled along the span, and the simulation was continued. 

This new mesh consisted of 500 x 100 x 53 points, and the initial condition corresponded to the steady 
solution from the previous simulation, replicated along the span, u contours were once again monitored, and 
the time series for this setup is shown in Fig. 3 

Similar to the previous setup, instabilities in the flow start to appear immediately, and after around 3 
FTTs, the flow starts to undergo numerical transition. The recycling/rescaling feed back process continues 
to increase the intensity of perturbations, which unlike in the narrower domain, become self-sustaining. Even 
after 15 FTTs, the perturbations are sustained, indicating that the flow has taken up a turbulent nature. 

To evaluate the response to even wider spans, the domain was once again doubled, resulting in a mesh 
with 500 x 100 x 105 points. The resulting L z = 240 w.u. The domain was initialized in a similar manner, 
using the statistically steady data from the L z = 122 w.u. simulation. The resulting mean velocity profiles 
at x = 12 Jr location are shown in Fig. 4. The relaminarized profile is indeed thinner than the two turbulent 
profiles, which have higher momentum content at lower y/5 locations. 

Having established the turbulent nature qualitatively, statistical averages of the velocity profile (Fig. 5) 
and Reynolds stresses (Fig. 6) were compared with DNS results of Spalart. These were computed at the 
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Figure 2. Time series of u (x z ,t) contours for a period of 15 flow-through times (FTTs) showing the apparent transition 
from laminar to a numerically turbulent flow, and subsequent relaminarization, for L z = 62 w.u. 


12£r location, where the simulated Re# = 500. Data was averaged along the spanwise direction, as well as 
over a number of flow through time cycles. 

To get an estimate on where the current simulation falls in comparison with the reference, spectral DNS 
mean velocity profiles at three different Re# locations (300, 670, 1,410) are plotted in Fig. 5 (a), along with 
the results from the implicit LES simulations. All three widths show perfect agreement within the linear 
viscous sublayer, until y + ~ 7. For the domain with L z = 62 w.u., the laminar profile is distinctly devoid 
of buffer and log layers thereafter, in addition to being offset in terms of the magnitude of the freestream. 
The turbulent profiles with L z = 122 and 240 w.u. wide domains were averaged over 20 FTTs, and more 
reasonable agreements are evident. The buffer and log layers become apparent, though the freestream 
magnitude is still relatively large. Whereas the ideal freestream for Re# = 500 profile should lie between 
the Re# = 300 and 670 curves, that for L z = 122 w.u. falls between Reg = 670 and 1,410 curves, and the 
L z = 240 w.u. freestream matches Reo = 670. These differences are attributed to the boundary condition 
set at the upper wall of the domain. By extrapolating the state variables, the condition at the upper wall 
is similar to that of slip-wall or symmetry boundary, and comparisons of the mean profiles do demonstrate 
that (not shown). Such a boundary condition causes the freestream to accelerate, resulting in higher mean 
velocities. 

Slight differences between spectral DNS and implicit LES are also evident within the buffer layer. This 
region is magnified in Fig. 5 (b). u + for both L z = 122 and 240 w.u. are marginally higher than the reference 
flat plate profiles of Spalart. This, again, is attributed to the symmetry-like boundary condition at the upper 
wall, which essentially makes the current results similar to that of a channel flow. For comparison, channel 
results of Moser et al. are also shown in the figure, and there is better agreement of the computed buffer 
layer with them. 
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(a) Initial 


(e) 6 FTTs 


(i) 12 FTTs 



(d) 4 FTTs 


(h) 10 FTTs 


(1) 15 FTTs 






Figure 3. Time series of u contours for a period of 15 FTTs showing the apparent transition from laminar to a 
numerically turbulent flow, for L z = 122 w.u. 


Reynolds stresses comparisons are shown in Fig. 6. <uu > , <vv>, <ww> and <uv> averages for 
L z = 122 and 240 w.u. are compared with reference DNS data. Since the L z = 62 w.u. case relaminarized, 
that is not shown here. For L z = 122 w.u., deviations from the references values are more pronounced, 
both, within the buffer layer and away from the wall (log layer). <uu> and <vv> fare poorly throughout 
the boundary layer, and it is concluded that the generated turbulence is non-physical for this span width. 
There is improved agreement for L z = 240 w.u., but since spectral DNS for Re# = 500 is unavailable, precise 
conclusions could not be drawn. Nevertheless, statistics within the viscous sublayer and the log layer show 
significant improvements, but there are differences evident within the buffer layer, especially for <uu> and 
<vv>. 

III.B. Re# = 1,400 

Since spectral DNS data for Re# = 1,410 is available, a more definitive comparison could be made for this 
case. To simulate Re# = 1,400 at a location 12 £r downstream of the inlet, Re# at the inlet was set at 1,048, 
based on Prandtl’s one-seventh power law approximation. 

Here again, an initial mesh of 500 x 100 x 27 points was considered. Based on the flow conditions, mesh 
resolution at the wall was computed to be 19.4 w.u. and 6 w.u. in the streamwise and spanwise directions, 
respectively. In the wall normal direction, the spacing was 0.65 w.u. at the wall, and was again stretched 
using a hyperbolic-tangent function. The overall span width was L z = 162 w.u. This is comparable to the 
resolution using by Spalart for his spectral DNS. 

At this width, transition to numerical turbulence occurred with the initial grid itself after approximately 
18 FTTs. Similar to the lower Re# case, this initial domain was doubled twice along the span, resulting in 
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Figure 4. Mean velocity profile for different span widths, at Re# = 500. 




(a) Viscous sub-layer, buffer layer and log-law region (b) Buffer layer, magnified 

Figure 5. Mean velocity profile in wall units for different span widths, at Re# = 500, compared with DNS results for 
fiat plate of Spalart, and channel flow of Moser et al. 
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Figure 6. Reynolds stresses for different span widths, normalized by friction velocity, at Reg = 500, compared with 
DNS results of Spalart 


domains that had 318 and 630 w.u. spans, respectively. The mean velocity profiles at the 12 £r location are 
shown in Fig. 7. That the boundary layer is thicker with increasing span widths is evident, with the mean 
velocities tending towards an asymptotic profile. 

Statistical averages of the mean velocity (Fig. 8) and Reynolds stresses (Fig. 9) at the 12c)r location are 
compared with reference DNS. Data was averaged over 20 FTTs for each setup, and statistical stationarity 
of the solution was established. 

In Fig. 8 (a), all three spans show evidence of the viscous sublayer, buffer layer and the log layer, unlike 
the previous case. There is, once again, very good agreement with spectral DNS for the viscous sublayer, until 

~ 10. The freestreams for all three cases are over-predicted, progressively dropping in magnitudes as the 
domains are widened. Deviations within the buffer layer are clearly evident here. This region is magnified in 
Fig. 8 (b), where spectral DNS of channel flow of Moser et al. are also included for clarity. The agreement 
of the current implicit LES with channel data can be clearly seen, confirming that the upper-wall boundary 
condition indeed influences the buffer zone. The acceleration of the freestream is also attributed to this 
symmetry-like boundary condition. 

For an ideal boundary layer setup, the top wall must accommodate the displacement in the freestream 
due to the presence of the boundary layer. A simple extrapolation does not sufficiently account for that 
displacement thickness. A vertical velocity profile, similar to that used by Lund et al., 18 determined from the 
integrated continuity relation, v(x) = U e dS* /dx + (5* — h)dU e /dx, should accommodate for that displacement 
thickness. In this study, since these are auxiliary simulations to channel-like configurations of separated flows, 
that boundary condition was not investigated. 

Comparisons of Reynolds stresses with reference DNS are shown in Fig. 9. Overall agreement significantly 
improves with increasing spanwise widths, with L z = 630 w.u. showing excellent agreement for <uu>, <ww> 
and <uv> profiles throughout the boundary layer. There is an over-prediction of the peak <uu> stress, 
which lies within the buffer zone, attributed once again to the upper-wall boundary condition. The influence 
of the higher peak <uu> is evident throughout the entire boundary layer. 
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Figure 7. Mean velocity profile for different span widths, at Re# = 1400. 


IV. Conclusions 

Direct simulations of zero-pressure gradient, turbulent boundary layers for two separated Re^ values 
(500 and 1,400) were conducted, and compared with spectral DNS data. The high-order compact-difference 
scheme and filter, with the standard recycling/rescaling method of turbulence generation, used meshes with 
similar resolutions as DNS, and the impact of varying span widths was investigated. 

A minimum span width was required for the mere initiation and continuation of Kelvin-Helmholtz-like 
numerical instabilities. Narrower domains (< 100 w.u.) failed to sustain the perturbations, resulting in 
relaminarization. Wider spans did sustain the instabilities. Comparisons of mean velocity profiles and 
Reynolds stresses showed that a span of 630 w.u. was required for the turbulent statistics to match reference 
DNS. The upper- wall boundary condition for this setup generated marginal deviations in the mean velocity 
and Reynolds stress profiles, particularly in the buffer zone. A symmetry-like upper-wall resulted in better 
agreement with channel data within the buffer zone, and a marginally faster freestream. 
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Figure 9. Reynolds stresses for different span widths, normalized by friction velocity, at Re# = 1,400, compared with 
DNS results of Spalart 
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